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The electrostatic continuum solvent model developed by Fattebert and Gygi is combined with a 
first-principles formulation of the cavitation energy based on a natural quantum-mechanical defini- 
tion for the surface of a solute. Despite its simplicity, the cavitation contribution calculated by this 
approach is found to be in remarkable agreement with that obtained by more complex algorithms 
relying on a large set of parameters. Our model allows for very efficient Car-Parrinello simulations 
of finite or extended systems in solution, and demonstrates a level of accuracy as good as that of 
established quantum-chemistry continuum solvent methods. We apply this approach to the study of 
tetracyanoethylene dimers in dichloromethane, providing valuable structural and dynamical insights 
on the dimerization phenomenon. 



I. INTRODUCTION 

The importance of electronic structure calculations in 
solution is self-evident: chemistry in nature and in the 
laboratory often takes place in water or other solvents, 
or at a solid-solvent interface. This is true for all of 
biochemistry, for most of organic, inorganic, and ana- 
lytical chemistry, and for a vast part of materials and 
surface sciences. The natural solution to this problem 
is to explicitly include the solvent molecules in the sys- 
tem, either as one or several solvation shells or as a 
bulk medium that fills the simulation box in periodic 
boundary conditions. Such approach rapidly increases 



the expense of the calculation and is not always afford- 
able. The reasons are twofold: the cost of an electronic- 
structure calculation scales as the cube of the number 
of atoms considered, at fixed density. Also, one needs 
to ensure that the solvent is treated appropriately as 
a liquid medium, using e.g. extensive Monte Carlo or 
molecular dynamics simulations. Given the large ratio 
between the number of degrees of freedom in the solvent 
vs. the solute, the statistical accuracy needed makes 
most of these approaches prohibitively expensive. The 
use of hybrid quantum- mechanical/molecular- mechanics 
(QM/MM) techniques^"— in which the solvent atoms 
are represented with point (or Gaussian) charges and 



classical potentials, can sensibly alleviate the cost of the 
computations, but does not remove the requirement of 
long dynamical trajectories of the combined quantum 
and classical fragments to simulate the liquid state of 
the solvent and to extract thermodynamical averages. 

Alternative to these explicit approaches, a description 
of the solvent as a continuum dielectric medium sur- 
rounding a quantum-mechanical solute has long been 
established, and has proved efficient and accurate in a 
diversity of cases In continuum schemes the dielec- 
tric fills the space outside a cavity where the solute is 
confined; the shape of this cavity, considered as a sin- 
gle sphere^ or ellipsoid in early implementations, has 
evolved to more realistic molecular shapes such as those 
defined by interlocking spheres centered on the atoms or 
by isosurfaces of the electron density.'''^ In the context of 
continuum models the interaction between the dielectric 
medium and the charge distribution of the solute provides 
the electrostatic part of the solvation free energy, AGei, 
which is the dominant contribution for polar and charged 
solutes. Solvation effects beyond electrostatic screening, 
conventionally partitioned in cavitation, dispersion, and 
repulsion, ^« are also important and will be discussed in 
the context of our model in Section II. In principle, the 
application of continuum models demands that no strong 
specific interactions are present between the solvent and 
the solute molecules, although the solvent can always be 
reintroduced explicitly as an "environmental" skin for the 
first solvation shells. 
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Inexpensiveness is not the single asset of continuum 
models against explicit solvent methods. Unless Monte 
Carlo or molecular dynamics techniques are used, it is 
unclear what orientation to choose for the molecules 
in an explicit solvent model, and even for a medium- 
sized solute there may be a large number of possible 
configurations with multiple local minima.'' More im- 
portantly, geometry relaxations will describe a solid or 
glassy phases for the solvent, with a mostly electronic di- 
electric screening that may differ substantially from its 
static limit. This is particularly true for water, where the 
static permittivity eq of the liquid is larger by a factor 
of twenty than its electronic too contribution. When ge- 
ometry optimizations including many solvent molecules 
are performed, changes in the solute — e.g. the hydration 
energy — remain "buried" or hidden by the large contribu- 
tions coming from the energy of the solvent. To extract 
meaningful information, Monte Carlo or molecular dy- 
namics simulations with accurate thermalizations and av- 
eraging times are necessary. Still, it is far from clear that 
even first-principles molecular dynamics treatments of a 
solvent would provide the accuracy needed to reproduce 
static screening as a function of temperature (as an exam- 
ple, the dielectric constant of water varies between 87.8 
at °C and 55.8 at 100 °C). Room temperature is well 
below the Debye temperature of many solvents, and thus 
the effect of quantum, Bose-Einstein statistics can be 
very important. In fact, recent first-principles molecular 
dynamics studies of water point to the fact that a combi- 



nation of inaccuracies in the quantum-mechanical models 
(such as density-functional theory in generalized-gradient 
approximations) and the use of Boltzmann statistics pro- 
duce an overstructured description of waXe'^~—, with 
apparent freezing roughly a hundred degrees above the 
experimental point. Last, the relaxation times needed 
to extract thermodynamical data from a solvated system 
can be exceedingly longfi^ compounding many of the is- 
sues highlighted here (dynamical, as opposed to static 
screening, would require to take into account the sol- 
vent relaxation times, either explicitly or via a frequency- 
dependent dielectric model, but such a framework goes 
beyond the scope of this paper). Continuum solvent 
methods are free from these issues, and for this reason 
alone they may be the first choice even when computa- 
tional resources are not the main constraint. 

The presence of a polarizable dielectric will induce a 
charge redistribution in the solute, which in turn will af- 
fect the polarization of the medium. In the self-consistent 
reaction field approach (SCRF) the dielectric medium 
and the electronic density respond to the electrostatic 
field of each other in a self-consistent fashion.* Over the 
past twenty five years a number of developments stem- 
ming from the SCRF approach have been proposed and 
further elaboratedii^"— Among these, the Polarizable- 
Continuum Model (PCM) of Tomasi et ali-^i^?i^? and the 
Conductorlike Solvation Model (COSMO) of Klamt and 
Schiiiirmannii are probably the most- widely used choices 
in quantum chemistry applications. In both cases the 
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dielectric constant e is taken to be 1 inside the cavity, 
and a fixed value outside (equal to the dielectric con- 
stant of the solvent for PCM, or infinite for the case of 
COSMO). The electrostatic problem is then formulated 
in terms of apparent surface charges (ASC) distributed 
on the solute-solvent interface. For first-principles molec- 
ular dynamics applications, the discontinuity of e at the 
interface needs to be removed to calculate accurately the 
analytic derivatives of the potential with respect to the 
ionic positions. This may be accomplished with the use of 
a smoothly varying dielectric potential that restores well- 
behaved analytic gradientsiSi Still, Born-Oppenheimer 
ab-initio molecular dynamics in localized basis sets are 
demanding enough that they have yet to be combined, 
to the best of our knowledge, with the ASC approach for 
realistic simulations of medium or large system. 

On the other hand, first-principles implementations of 
the continuum solvent model within the Car-Parrinello 
frameworli24 have been devisedf2k~— even though dy- 
namical studies have been reported, to the best of our 
knowledge, in only few cases^SiSl In this paper, we intro- 
duce a first-principles and conceptually simple approach 
to the calculation of cavitation energies based on the def- 
inition of a quantum surface for the solvent^ We com- 
bine this scheme with the electrostatic solvation model of 
Fattebert and Gygi(2fi*29^ g^^^ gj-^j level of accuracy at 
least as good as that of established quantum-chemistry 
treatments. The model requires no adjustable parame- 
ters other than a universal definition of the cavity (prac- 



tically depending on one parameter), and the dielectric 
constant and the surface tension of the solvent. This 
combined model is well suited for first-principles molec- 
ular dynamics calculations of large finite and extended 
systems, using e.g. efficient plane-wave Car-Parrinello 
implementations. In the following sections we describe 
the method and examine its performance in compari- 
son with experiments and with the well-established PCM 
approach. Finally, given that cavitation contributions 
can be particularly important in dimerization processes 
(where the fusion of two cavities into one provides an 
additional stabilizing energy), we employ our method to 
study the association of the tetracyanoethylene (TCNE) 
anion in solution'^'' by means of static and dynamical sim- 
ulations, highlighting the role of the cavitation term in 
the dimerization. 

II. THE MODEL AND ITS CONTEXT 
A. Preliminary details 

Our continuum solvation model has been implemented 
in the public domain Car-Parrinello parallel code in- 
cluded in the Quantum-ESPRESSO package;^ based on 
density- functional theory (DFT), periodic-boundary con- 
ditions, plane-wave basis sets, and pseudopotentials to 
represent the ion-electron interactions. All calculations 
reported in this work, unless otherwise noted, have been 
performed using Vanderbilt ultrasoft pseudopotentials^ 
with the Kohn-Sham orbitals and charge density ex- 



4 

panded in plane waves up to a kinetic energy cutoff of 
25 and 200 Ry respectively. In the Appendix we review 
the formalism used to calculate energies and forces in pe- 
riodic boundary conditions in the context of our imple- 
mentation. Further details can be found in reference—. 

We adopt the definition introduced by Ben-Naim for 
the solvation free energy, '^'^ in which AGsoi corresponds 
to the process of transferring the solute molecule from 
a fixed position in the gas phase to a fixed position 
in the solution at constant temperature, pressure, and 
chemical composition. For calculation purposes and es- 
pecially in the case of the continuum dielectric model, 
AGsoi can be regarded as the sum of several compo- 
nents, of which the electrostatic, the cavitation, and the 
dispersion-repulsion contributions are the most relevant 
(AGsoi = AGe( -1- AGcav + AGMs-rep)<^ Nouc of thcsc, 
however, can be directly obtained through experiment, 
the sum of all of them, AG^oi, being the only measur- 
able quantity. In our model, AGei and AGcav are con- 
sidered explicitly, while AGdis~rep, less relevant for the 
systems considered here, is largely seized by virtue of the 
parametrization, as part of the electrostatic term. The 
dispersion-repulsion energy may be important in the case 
of hydrophobic and aromatic species, but its explicit cal- 
culation is beyond the aim of the present work — in par- 
ticular, the implementation of the technique proposed by 
Floris, Tomasi and Pascual Ahuiri2L2& would be straight- 
forward in our model. 



B. Electrostatic solvation energy 

The electrostatic interaction between the dielectric and 
the solute is calculated as proposed by Fattebert and 
Gygii^S*^ In the following we provide an outline of the 
model. 

The Kohn-Shani energy functional^^ of a system of ions 
and electrons can be written as 

E[p\=T[p\ + J v{v)p{v)dr + E,, + ^ J p{r)cb[p]dr (1) 

where the terms on the right hand side correspond to 
the kinetic energy of the electrons, the interaction energy 
with the ionic potential, the exchange-correlation energy, 
and the electrostatic energy E^s respectively. In the stan- 
dard energy functional, the electrostatic potential is 
the solution to the Poisson equation in vacuum, 

V2(/) = -47rp . (2) 

In the presence of a dielectric continuum with a permit- 
tivity e[p], the Poisson equation becomes 

V • (e[p]V0) = -iirp . (3) 

By inserting the charge density obtained from Eq. (3) 
into the expression for the electrostatic energy, and inte- 
grating by parts, we obtain: 

= ^1 e[p]{Vcj,[p]fdr. (4) 

While Eq. (2) can be efficiently solved in reciprocal space 
with the use of fast Fourier transforms, for arbitrary e[p] 
the Poisson equation (3) must be solved with an alter- 
native numerical scheme. In the present case, it is dis- 
cretized on a real space grid, and solved iteratively using 
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a multigrid technique.— The functional derivative of E^s 
with respect to p yields (p and an additional term V^, 
originating in the dependence of the dielectric function 
on the charge density: 

^(r)=0(r) + y,(r), (5) 

^e(r) = -^(V</)(r)f^(r). (6) 

The self-consistent Kohn-Sham potential is constructed 
summing and the electrostatic potential 0, to which 
contributions from the exchange-correlation, and the lo- 
cal and non-local terms in the pseudopotentials are also 
added (see Appendix). The dielectric medium and the 
electronic density then respond self-consistently to each 
other through the dependence of e on p and viceversa. 

As already mentioned in the introduction, in Gaussian- 
basis sets implementations of the continuum model e is a 
binary function with a discontinuity at the cavity surface. 
The accurate representation of such a function would re- 
quire unrealistic high kinetic energy cutoffs for the plane 
wave basis and expensive real space grids. The use of 
smoothly varying dielectric functions instead eases the 
numerical load and avoids discontinuities in the forces, 
essential to proper energy conservation during molecular 
dynamics simulations. Also, a smooth decay of the per- 
mittivity in the proximity of the solute-solvent boundary 
may even be considered a more physical representation 
than a sharp discontinuity. In our implementation the di- 
electric medium is defined using two parameters po and 



/3: 

This function asymptotically approaches Coo (the permit- 
tivity of the bulk solvent) in regions of space where the 
electron density is low, and 1 in those regions where it is 
high. The parameter po is the density threshold deter- 
mining the cavity size, whereas f3 modulates the smooth- 
ness of the transition from eoo to 1. 

C. Cavitation energy 

The cavitation energy AGcav is defined as the work 
involved in creating the appropriate cavity inside the 
solution in the absence of solute-solvent interactions)^ 
Different approaches have been introduced to compute 
AGcav', nevertheless it is unclear which one is the most 
accurate given the unavailability of experimental values 
to compare. Formulations based on the scaled particle 
theorji^Sdi have been originally proposed by Pierotti-- 
and further developed in several studiesi^"— Although 
these approaches are derived from a rigorous statistical 
mechanics standpoint, eventually the use of a set of fit- 
ted parameters is needed to represent an effective radius 
for the solvent and for the spheres centered on the so- 
lute atoms. For nonspherical cavities, one of the most 
used approximations is the so-called Picrotti-Claverie 
formulatSi^ 

AGcav = ^ ^^„2 Gcav{Rk)- (8) 
k=l 
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Eq. (8) describes the cavity as the volume occupied by 
N interlocked spheres centered on the atoms; At is the 
area of atom k exposed to the solvent, Rk is its van-der- 
Waals radius, and Gcav{Rk) is the cavitation free energy 
associated to the creation of a spherical cavity of radius 
Rk according to Pierotti4^ 

Efforts have also been made to describe AGcav as a 
function of the macroscopic surface tension of the sol- 
vent 7ii&~iSi The suggestion of Uhlig^ of expressing the 
work involved in producing the cavity as the product be- 
tween 7 and the area of a sphere, AGcav ~ 47ri?^7, has 
been extended to account for the curvature of the solute- 
solvent interface, according to the theory of Tolman for 
the surface tension of a dropleti^ The validity of simpli- 
fied expressions of the kind 

AGcav =PV + AnR^ (^1 - (9) 

has been investigated by different authors^^-^'^ by means 
of Monte Carlo simulations with classical potentials. In 
Eq. (9), 7 is an effective surface tension for the interface, 
R is the radius of the cavity, and (5 is a coefficient that 
would correspond to the Tolman length in the case of 
a macroscopic surface. Studies from both Florist and 
Chandlcr^'^ groups have shown that 7 is essentially in- 
distinguishable from the macroscopic surface tension of 
the solvent, 7. Their simulations have assigned to (5 a 
value of 0.0 in TIP4P water and of the order of -0. Sa- 
in the case of different Lennard- Jones fluids (cr being the 
Lennard- Jones radius) suggesting that the curvature 
correction can in practice be ignored for cavities with 



radii above only a few Angstroms. 

In view of these results, we have chosen to estimate 
the cavitation energy as the product between the surface 
tension and the area of the cavity, 



(10) 



where S{pq) is the surface of the same cavity employed 
in the electrostatic part of the solvation energy and is de- 
fined by an isosurface of the charge density. As observed 
by Floris et al.jSi there is always a surface in between 
the internal and the solvent accessible surfaces such that 
the correction factor (1 — ) reduces to 1, entailing a 
linear dependence between AGcav and the cavity area. 
We rely on the parametrization of the density threshold 
Pq to obtain an appropriate surface. 

The area of this cavity can be easily and accurately 
calculated by integration in a real-space grid, as the vol- 
ume of a thin film delimited between two charge density 
isosurfaces, divided by the thickness of this film. This 
idea has been originally proposed by Cococcioni et alMi 
to define a "quantum surface" in the context of extended 
electronic-enthalpy functionals : 

|Vp(r)| 



S{po) = J dr{^^^_. [p(r)] [p(r)]} 



(11) 



The finite-differences parameter A determines the sep- 
aration between two adjacent isosurfaces, one external 
and one internal, corresponding to density thresholds 
Pa — A/2 and pa + A/2 respectively. The spatial dis- 
tance between these two cavities — or the thickness of the 
film — is given at any point in space by the ratio A/|Vp|. 



The (smoothed) step function is zero in regions of low 
electron density and approaches 1 otherwise, and it has 
been defined consistently with the dielectric function of 
Eq. (7): 



{p{v)lpa?P ~ 1 



1 



(12) 



.(p(r)/po)2^ + l 
Note that the volume of the cavity is simply the integral 
of 1? on all space: 



Vc{po) 



dri?pjp(r)]. 



(13) 



The functional derivative of lS.Gcav — lS{p) with respect 
to the density gives then the additional contribution to 
the Kohn-Sham potential, 

= Jx[,?^„_.[p(r)]-^,„^.[p(r) 



V^V^ dip{v)d.jp{v)didjp{Y) _ ^ 

i j 



dfp{r) 



(14) 



|Vp(r)|3 ^ |Vp(r)| 

where the indices i and j run over the x, y, z coordinates, 
and di indicates a partial derivative with respect to the 
position. 

The exact value of the discretization A is not impor- 
tant, as long as it is chosen within certain reasonable 
limits — a very low value would introduce numerical noise, 
while a too large one would render an inaccurate mea- 
sure of the surface. The freedom in the choice of A is 
illustrated in Fig. where the dependence of S on this 
parameter is examined for a water molecule at various 
thresholds. For po equal or above 0.00048 e, the calcu- 
lation of the cavity area is fairly converged for any value 
of A within the range displayed in the figure. We have 



adopted a value of A=0.0002 e in our simulations. It is 
worth noting, on the other hand, that the dependence of 
the surface on the density threshold po is only moderate, 
reflecting the fact that at the "molecular boundary" , the 
electron density decays significantly on a short distance. 
This behavior is portrayed in Fig. ^ where it can be seen 
that for a given A, the calculated surfaces change in only 
about 25% when po is increased three times. AGcav is in 
fact much less sensitive to the electron density threshold 
than AGei- 

III. RESULTS AND DISCUSSION 
A. Solvation energies in water 

The only adjustable parameters in our solvation model 
are po and (3, which determine the shape of the cavity 
according to Eqs. (7) and (12). Other parameters enter- 
ing the model, namely the static dielectric constant and 
the surface tension of the solvent, are physical constants 
taken from experiments. We have actually kept po as 
the single degree of freedom to fit the solvation energy, 
while fixing the value of /3 to 1.3 as in reference^. This 
choice of P provides a smooth, numerically convenient 
transition for the step function, still ensuring that the 
lower and upper limits oi e(p{r)) and •&{p{r)) are reached 
reasonably fast. The parameter po was obtained from a 
linear least squares fit to the hydration energies of three 
solutes: amide, nitrate, and methylammonium (a polar 
molecule, and two ions of opposite sign). The resulting 
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value, Po =0.00078, was employed thereafter in all the 
simulations. This can be regarded as a rather universal 
choice for po and /3; reparametrizations for different sol- 
vents could be considered (if enough experimental data 
were available) probably gaining some marginal accuracy 
at the expense of generality. 

Table I shows the solvation and cavitation energies in 
water calculated for a number of neutral species, along 
with their experimental values^"— A quite remarkable 
agreement with experiments is found. We compare the 
data with PCM results obtained at the DFT-PBE/6- 
311G(d,p) level (or DFT-PBE/3-21G** for the case of 
Ag+) using the Gaussian 03 packaged Also significant 
is the accord between the cavitation energies computed 
with the two methods — with the caveat that in Gaussian- 
PCM AGcav is based on the Pierotti-Claverie formula 
(see Eq. (8)) which requires a lengthy list of parame- 
ters including all van-der-Waals radii. Similar agreement 
between the values of AGcav coming our approach and 
PCM is found among charged solutes as shown in Ta- 
ble II. The level of accuracy in AGsoi is in this case as 
good as for the neutral solutes, if viewed in relative terms 
(we point out that, regarding the experimental values of 
AGsoi reported for ions, discrepancies between sources 
up to a few kcal/mol are common). 

The solvation energies of the ionic solutes showed 
in Table II were calculated including the Makov-Payne 
correctioufSi which takes into account how the gas phase 
energy of a charged system is affected by its periodic im- 



ages in supercell calculations: 

EGAS = EpBC+''^-'^-^ + 0[L-% (15) 

where Eg as and EpBC are the isolated and the supercell 
energies respectively, q is the charge of the system, Q 
its quadrupole moment, L the lattice parameter, and a 
the Madelung constant (we used a simple cubic lattice 
of charges, for which a— 2.8372^'^). As shown in Fig. [21 
for the nitrate anion, the dependence of the energy with 
respect to the inverse of the lattice parameter becomes 
virtually linear for L above 40 a.u., pointing out that the 
quadrupole term can be neglected in supercells of that 
size or larger. So, we applied the Makov-Payne correction 
to the 1 /L leading order to all the cations and anions in 
Table II, always checking for convergence with respect to 
1/L. The gas phase energies calculated in this way were 
subtracted from the correspondent energies in solution 
to obtain AG sol- Fig. |21 also shows that total energies 
in solution quickly converge with respect to the size of 
the supercell, thanks to the dielectric screening of the 
Coulombic interactions between periodic images. 

With the exception of CH4 , the solutes in Tables I and 
II are either polar or ionic compounds of relatively small 
size. Since the dispersion-repulsion energy is not explic- 
itly accounted for in our model, its accuracy for systems 
in which this contribution becomes dominant, such as 
highly hydrophobic or aromatic compounds, will be nec- 
essarily affected (this is already the case for methane). 
For the species listed in Tables I and II the dispersion- 
repulsion effect is captured to a large extent by our elec- 
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trostatic term. As the size of the solute increases and its 
polarity decreases, though, the non-electrostatic terms 
tend to monopolize the solvation energy, and a model 
lacking the dispersion-repulsion contribution will perform 
poorly. This limitation could be possibly overcome by a 
different parametrization specific to large nonpolar so- 
lutes, or of course by directly computing the dispersion 
and repulsion contributions. 

B. Molecular dynamics and total energy 
conservation 

As a computational test, we performed canonical 
molecular dynamics simulations of a tetramer of D2O 
molecules (heavy water), using our solvation-cavitation 
model to represent an aqueous solution. A time step of 
0.192 fs and an electronic mass of 400 a.u. were employed, 
and the system was thermalized at 350 K by applying the 
Nose-Hoover thermostat on the ions. Fig. O shows the 
initial configuration of the cluster, where the four D2O 
molecules are stabilized in a ring by four hydrogen bonds. 
In the upper part of Fig.^ the total energy is monitored 
throughout the run and compared with the potential en- 
ergy. The conservation of the total energy is as good 
as in the gas phase for the same simulation parameters, 
and is not affected by the dissociation of the bonds. The 
analysis was not pursued beyond 0.65 ps, when the wa- 
ter cluster dissolves into the medium and one of the D2O 
molecules evolves close to the border of the real space 
grid, affecting the Dirichlet boundary conditions. 



In the lower section of Fig.01 the intermolecular O - • -H 
distance is plotted for the four initial hydrogen bonds 
that keep the cluster bound. The solvent dissociates this 
structure early in the simulation, and before half a pi- 
cosecond only one hydrogen bond has survived. By the 
end of the run a dimer is what remains of the original 
tetramer. For comparison, long molecular dynamics (up 
to 10 ps) were carried on in the gas phase under identical 
conditions. In this case the cyclic cluster is stable for the 
full length of the simulation, showing that the disruption 
of the intermolecular bonds is indeed a consequence of 
the solvation effect. 

C. Dimerization of TCNE anions in solution 

Starting in the early 60 's, dimerization of charged and 
neutral organic tt radicals in solution and in the solid 
state was reported by several authors The discov- 
ery of this phenomenon prompted a vast amount of re- 
search which has continued up to the present daymS^~— 
Among the systems addressed, significant efforts have 
gone into the study of the tetracyanoethylene anion 
[TCNE]^' and its salts because of their central role in 
the understanding and development of molecular met- 
als. Recently, evidence has been presented showing that 
the dimerization of [TCNE]^' in the solid state involves 
two-electron four-center 7r*-7r* bonding arising from the 
interaction of the two singly occupied molecular orbitals 
(SOMOs) of the anions and leading to long (w 3.0 A) 
intermonomer C-C covalent bonds. ^^'^^ Data from UV- 
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vis and EPR spectroscopies suggested the same con- 
clusions are true in the solvated statei^SiSS In the gas 
phase, DFT and MPS calculations show that the dimer 
is only metastable, since the attractive covalent interac- 
tion between the anions is outweighed by the Coulom- 
bic repulsioniSSiSSiZS. In the solid state, in contrast, the 
positive counterions stabilize the array of like charges, 
allowing the 7t*-tt* bonding to occmM^ 

A solvent may play an analogous role in stabilizing the 
dimer, by favoring the concentration of charge in a single 
cavity. We used our solvation model to fully optimize 
the doubly charged TCNE dimeuS in dichloromethane, 
properly adapting the values of e and 7 (8.93 and 27.20 
mN/m respectively). We found a stable minimum at an 
equilibrium distance of 3.04 A, in close agreement with 
solid state geometries: X-ray data of different salts2^~ii& 
range from 2.83 to 3.09 A. The CN substituents devi- 
ate from the plane by 5° (see Fig. , consistently with 
the NC-C-C-CN dihedral angles observed in crystals, be- 
tween 3.6 and 6.5°. This deviation has been ascribed to 
the rehybridization of the sp^ carbon as the intradimer 
bond is formed, but its origin could be also tracked to 
the steric repulsion between the CN moieties facing each 
other. 

Fig. 13 (upper panel) shows the binding energy for the 
[TCNE]2^ in dichloromethane as a function of the sepa- 
ration between the [TCNE]^' fragments. At every point, 
all coordinates were relaxed while freezing the intradimer 
distance. The curve presents a steep minimum at 3.04 A, 



with a barrier to dissociation of nearly 4 kcal/mol. The 
grouping of two monomers inside a single cavity, of an 
area smaller than the one corresponding to two separate 
cavities containing one monomer each, is energetically 
favored by the surface tension of the solvent. Thus, if 
the contribution of the cavitation energy to the solvation 
is not considered, the binding results weaker, as seen in 
Fig. El The surface of the cavity, plotted in the lower 
panel, increases gradually as the monomers are pulled 
apart, until the solvation cavity splits in two at around 5 
A. (this is a case in which different (3 in the parametriza- 
tion could account for the distinctive ability of solvents 
to penetrate narrow spaces). Beyond this point the total 
surface remains constant as each [TCNE]^' unit occupies 
a separate cavity, and the two curves in the top panel 
merge. The ground state of the system is a singlet for 
distances up to 4.0 A, whereas at larger separations the 
spins of the fragments are no longer paired, conforming 
to a triplet state. 

A value of -1.1 kcal/mol is obtained for the binding 
energy between the monomers. Such a value is underes- 
timated with respect to the experimental dimerization 
enthalpy, AHd, reported in the range of -6.9 - -9.8 
kcal/mol in dichloromcthaneiiSii The disagreement can be 
partially attributed to the inability of DFT to fully ac- 
count for the correlation energy involved in the 7r*-7r* 
bond, and also, to some extent, to the effect of the coun- 
terions present in the solution, which differentially stabi- 
lize [TCNE]^^ compared to two [TCNE]"- anions. This 
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effect has been advocated in a recent study*^ of the in- 
teraction of two [TCNE]^' fragments in tetrahydrofuran 
(e=7.58) using PCM at the MP2 level, to explain why the 
dimer was found metastable by 9.7 kcal/mol with respect 
to the isolated monomers — the experimental estimate for 
AHd being -8 kcal/mol in 2-methyl-tetrahydrofuranifi^ 
The binding energy curve presented in that work exhib- 
ited a broad minimum extending from 3.1 to 3.7 A, a 
separation range substantially larger than the one ob- 
served in the solid state. Our own PCM calculations in 
dichloromcthanc, using PBE in combination with the 6- 
311+G(d,p) Gaussian basis set, yield a metastable dimer 
with an interaction energy of 3.2 kcal/mol and an equi- 
librium distance of 3.00 A. 

Temperature dependence investigations in solution in- 
dicate that the dissociated [TCNE]^' anions arc the pre- 
dominant species at ambient conditions, and that the 
concentration of the dimer rapidly grows as the temper- 
ature goes down^S*^ Car-Parrinello molecular dynamics 
simulations of the [TCNEjj" dimer were performed in 
dichloromethane at 250 K, with the temperature con- 
trolled by the Nose-Hoover thermostat on the ions. A 
time step of 0.288 fs and an electronic mass of 400 a.u. 
were used. In Fig. 13 we monitor the evolution of two 
structural parameters which serve as descriptors of the 
[TCNE]- - [TCNE]~- bonding. The intradimer separa- 
tion, departing from a value of 3.9 A corresponding to an 
initially elongated dimer, drops to nearly 2.7 A and then 
describes large oscillations in the order of 1 A around the 



equilibrium distance. The second parameter, correspond- 
ing to the C=C- • -C^C dihedral angle formed by the two 
[TCNE]^' anions, provides a measure of the alignment 
between the monomers: if this angle is zero the anions 
lay parallel. Fig. [7| shows that this is not the case most 
of the time. Rapid oscillations of an average amplitude 
of 6° take place around the equilibrium angle. During 
most of the second part of the run the oscillations are 
not necessarily centered around zero, which is indicative 
of the relatively lax nature of the bond. 

The length of the simulation is enough to reveal some 
distinctive features of the frequency spectrum of the sys- 
tem in the IR region. The continuous line in Fig.|Hlshows 
the Fourier transform of the velocity- velocity correlation 
functions corresponding to two pairs of atoms in the 
dimcr. The first pair consists of the two carbon atoms in- 
volved in the C=C bond. The autocorrelation function of 
the relative velocity between these two centers originates 
an intense peak corresponding to the C=C stretching at 
1250 cm^^. The same mode resolved in the case of the 
monomer (dashed line) shows up at 1310 cm^^. In the 
solid state, experimental C=C stretching frequencies of 
1364 and 1421 cm~^ have been reported for the dimcr and 
the monomer respectively^ Such discrepancies between 
our results and the experimental numbers are expected, 
given the distinct conditions in the solid and liquid envi- 
ronments, the difference in the temperatures at which the 
spectroscopic and the computational data were collected, 
the use of DFT, and the slight downshift in ionic frequen- 
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cies in Car-Parrinello dynamics.^^ However, we note that 
the shift of 60 cm^^ in going from the monomer to the 
dimer is nicely reproduced by our simulations. 

In an attempt to characterize the frequency of the 
intradimer 7r*-7r* bonding, we have also analyzed the 
relative- velocity autocorrelation function for the two car- 
bon atoms forming the bond, one atom pertaining to 
each monomer. The frequency spectrum of this func- 
tion yields the four groups of signals appearing below 
600 cm^^ in Fig. |H1 the assignment of which is less evi- 
dent than in the case of the C=C stretching. Although 
we are unable to unambiguously identify all these fre- 
quencies, Fourier transform analysis of the autocorrela- 
tion function for the velocity of the center of mass of 
the two fragments (data not shown) points to the lowest 
frequency emerging in the spectrum, at 65 cm^^, as the 
one related to the intradimer vibration. To the best of 
our knowledge, no experimental data is available for this 
mode. Interestingly enough, though, the aforementioned 
theoretical study based on PCM and MP2,^^ predicted 
an inter- fragment vibrational frequency of 60 cm^^ by 
solving the one-dimensional Scrodinger equation on the 
potential energy surface calculated for the interaction be- 
tween the [TCNE]~- anions. 

IV. FINAL REMARKS 

The electrostatic-cavitation model described in this 
work enables Car-Parrinello molecular dynamics simula- 
tions in a continuum solvent for large finite systems, and 
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shows a level of accuracy as good as that offered by state- 
of-the-art quantum chemistry solvation schemes. Addi- 
tionally, our model is suited for the treatment of periodic 
systems in solution, representing a powerful tool for the 
study of solid-liquid interfaces, solvated polymers, and 
in general extended systems in contact with a solution. 
Further improvements will be the subject of future work, 
especially the incorporation of the dispersion-repulsion 
effects, which become increasingly important with the 
size of the solute. The method of references — and — is 
an attractive choice, although other possible approaches 
derived from first-principles and employing a minimal 
number of parameters are also envisioned. 

Our cavitation energy, defined in a simple and physi- 
cal way, can be straightforwardly implemented in plane- 
waves or real space codes. Interestingly, such definition 
turned out to be in remarkable agreement with the val- 
ues provided by more complex algorithms reliant on large 
sets of parameters. 

The real time study of the pairing of [TCNE]^' consti- 
tutes the first dynamical ab-initio investigation of dimer- 
ization phenomena in solution, of which the formation 
of the [TCNEjj^ is just one example. The binding of 
charged radicals in solution is relevant to a broad field of 
research in organic and materials chemistry, and proper 
consideration of the cavitation contribution turns out to 
be a central ingredient for an accurate atomistic descrip- 
tion. 
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VI. APPENDIX 

We sumarize here the relevant steps to calculate en- 
ergies and forces in the framework of pseudopotential 
codes in periodic boundary conditions, highlighting the 
additional terms arising from the electrostatic embed- 
ding. Leaving aside the exchange-correlation energy and 
the non-local term of the pseudopotential, the electro- 
static problem in a system of pseudo-ions (nuclei plus 
core electrons) and valence electrons may be written^ 

^ = E §^ + E / Peir)vioc{r - Ri)dr 

The first term on the right in Eq. (16) accounts for re- 
pulsion between pseudo-ions, the second is the interac- 
tion between these ions and the valence electron density, 
and the third is the Coulombic integral between valence 
electrons. Let pj{r — Rj) be a Gaussian distribution of 
negative sign that integrates to the total charge of the 
pseudo-ion (note that the electronic charge is defined here 
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as positive). Adding and subtracting J2i Pi{'''—Ri) from 
Peir) in the third term we obtain 

E=\f f[Pe{T) + Y,Pi{r-Ri)][Pe{r') 



ij 



+ ^ / Peir)vUr - Ri)dr + E (1^) 

The first term on the right is the Hartree energy Eh of 
a pseudopotential code. Introducing the following defini- 
tions: 



it is possible to write the total energy as: 

E = Eh + Eps + Egr + Egeif (18) 
The pseudo-ions density p/ is defined as: 

where Rj determines the width of the Gaussian associ- 
ated with the site I. Under such definition Eseif and 
Esr can be evaluated analytically. In particular, Esei / is 
a constant not dependent on the atomic positions: 



El \ ^ ZlZj 

Esr = — — er/c 



Rij 



(20) 



(21) 



Rij K^/iWrm? 
The remaining pseudopotential term Eps is computed in 
reciprocal space, after constructing the pseudopotential 
^z^oc carrying both contributions from the local pseudopo- 
tential vioc and the smeared core charges potential u/. 

Eps=J2j Pe{r)vUr)dr (22) 



-^ps = X] j Pe{'r)['^ioc{r - Ri) +vi{r- Ri)]dr, 



with vi{r) 



^dr' 



■E 

KJ 



ZiZj 
Rij 



'^loci^) = vioc{r) + vi{r) = vioc{r) - / \ dr' 



\r — r'\ 



..oc(r)-^er/(^ 



(23) 



The ionic forces can be obtained from the energies 
above (plus the non-local pseudopotential term, which 
will be omitted for simplicity). The Hellmann-Feynman 
theorem — i.e. the stationariety of the total energy with 
respect to ijj — gives 

SE Sltpj) _ dE 



p _ _dE _ _dE 
^ ~ ~dRi ~ ~dR'i ~ ^ S\^j) dRi ^ dRi 



(24) 
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(real wavefunctions are assumed). Thus, 

dEff dEps dEsr dEseif 



dE 
' dRi 



dE 



dRi dRi dRi dRi 



whose derivative with respect to the atomic positions is: 



(25) \ j j[^T.P^^r-R,)[^-^^Y.P^{r'-R,) 



Note that the partial derivatives of the individual terms 
in the Hamiltonian do not correspond to the total deriva- 
tives. For example: 



dEH{Ri, R2, ...,Rn) , dE 
mi, ^ dRi 



Um ^^(-^1' -^^5 -"J -Rn + e) — Eh{Ri,R2, Rn — e) 



Eaei f does not depend on Rj and therefore does not con- 
tribute to the forces, whereas the derivative for Esr can 
be obtained analytically. 



The derivative of Epg results 



a^=i:y^e(r)^i^c^r 



(26) 



where the term dvi^^{r)/dRi is straightforward in the 
reciprocal space: 



vUr-Rj)=Y,VGe'''^e-'' 



GRi 



+2pe(r)(^^P7(r'-i?/) 



r — r 



jdrdr' 



= //(al7?"<^'-'">) 

X {pe{r) +Y,pi{r- Ri)j y^^drdr' 

Hence, if ptot('') = Pe{r)+J2i Pi{^—Ri)-: the contribution 
from Eh turns out to be 

where the term d'^jpi{r' — Ri)/dRi is obtained in 
Fourier space in the same fashion as in Eq. (27). The 

ratio / dr' ptot{'r)l\r — r'\ is the Hartree potential Vh, 
which can be computed in the reciprocal space from the 
expansion for ptot{r)- 



Ptotir) = pee''"', Vh = J2 ^ce' 



Gr 



Vlo, 



OR 



{r - Ri) =J2-iGiGe'^''e- 



iGRi 



(27) 



with vg the coefficients of the Fourier expansion for 

Finally, to obtain the contribution from Eh, the 
Hartree energy is recast as: 

1 



Eh = 



.ir)peir')+Ypi{r-Ri)J2piir'-Ri) 



+2p,{r)Y,Pi{r' -Ri) 



r — r 



rdrdr' 



-47r 



^Pce 



iGr 



(29) 



In the case of the continuum solvent implementation, 
Vh is replaced by ^^(r) according to Eq. (5) and (6) 
of the main text. The electrostatic contribution to the 
energy originated in the dielectric medium is computed 



as 



E,_. 



j Pe{r)(f){r)dr 



(30) 



where ^(r) is the electrostatic potential obtained using 
the multigrid in Eq. (3). The Hartree term Eh is thus 
replaced by Egg in the calculation of the total energy. 
The cavitation energy is accounted properly in the to- 
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tal energy by adding 'yS, which functional derivative 
(Eq. (14)) is included in the Kohn-Sham potential — as 
the Hellmann-Feynman theorem applies for that term. 
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TABLE I: Solvation and cavitation free energies (kcal/mol) 
for neutral solutes in water, calculated with this model and 
with PCM as implemented in Gaussian 03. 



Expt!^ — This model PCM This model PCM 



iJ-2 ^ -' 


-fi S 
KJ.O 


8 4 




^ 7 


0. 1 


NH3 


-4.3 


-3.2 


-1.6 


6.6 


6.6 


CH4 


2.0 


5.4 


6.9 


7.5 


10.0 


CH3OH 


-5.1 


-3.6 


-0.8 


9.0 


9.6 


CH3COCH3 


-3.9 


-1.7 


3.5 


13.7 


14.3 


HOCH2CH2OH 


-9.3 


-9.3 


-6.7 


13.0 


12.3 


CH3CONH2 


-9.7 


-10.5 


-4.6 


12.7 


12.8 


CH3CH2CO2H 


-6.5 


-6.0 


-2.4 


14.8 


14.6 


mean unsigned error 




1.5 


4.0 






max. unsigned error 




3.4 


7.4 
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TABLE 11: Solvation and cavitation free energies (kcal/mol) 
for ionic solutes in water, calculated with this model and with 
PCM as implemented in Gaussian 03. 











^Gcav 




Expt^-^ 


This model 


PCM 


This model PCM 


cr 


-75 


-66.9 


-72.6 


7.9 


5.8 


NO3- 


-65 


-57.8 


-62.6 


10.5 


9.7 


CN" 


-75 


-64.8 


-70.2 


8.4 


7.0 


CHCl2CO;7 


-66 


-74.7 


-53.5 


16.3 


15.7 


Ag+ 


-115 


-110.0 


-102.3 


5.7 


4.0 


CH3NH+ 


-73 


-81.0 


-65.1 


9.4 


10.2 


CH3C(0H)CH+ 


-64 


-70.6 


-55.2 


13.5 


14.4 


C5H5NH+ (pyridinium) 


-58 


-60.8 


-59.0 


15.0 


13.9 


mean unsigned error 




7.1 


6.6 






max. unsigned error 




9.2 


12.7 
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Figure Captions: 

Figure 1.: Cavity area of a water molecule as a function 
of A (thickness parameter used to evaluate the area, see 
text) for several values of the electronic density threshold 
Po- 

Figure 2.: Total energy of the NO3 anion as a function 
of the inverse of the lattice parameter, computed in vac- 
uum, in solution, and in vacuum with the Makov-Payne 
correction up to the leading order. 

Figure 3.: Cluster of D2O molecules used as start- 
ing configuration in the molecular dynamics simulations 
which results are reported in Fig. 4. 

Figure 4.: Total and potential energies (top) as a func- 
tion of time in a molecular dynamics simulation of a cyclic 
tetramer of heavy water in aqueous solution. The to- 
tal energy contains the contribution of the Nose- Hoover 
thermostat. The four curves starting at the bottom of 
the graph represent the evolution of the intermolecular 
O- • -H distance between the atoms initially involved in 
hydrogen bonds. 

Figure 5.: Optimized structure of a dimer of [TCNE]"' 
in dichloromethane, enclosed by an electronic density iso- 
surface at 0.00078 e delimiting the solvation cavity. Car- 
bon atoms in light gray and nitrogen atoms in dark. 

Figure 6.: Upper panel: binding energy of two 
[TCNE]"' anions in dichloromethane as a function of its 
separation, calculated with only the electrostatic contri- 
bution to the solvation energy, and with both the electro- 
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static and cavitation contributions. Lower panel: area of 
the solvation cavity as a function of the separation be- 
tween the [TCNE]"' anions. Above 5 A the cavity splits, 
and the plotted values correspond to the area of two cav- 
ities containing one [TCNE]"' each. 
Figure 7.: Time evolution of the intradimer separation 
(top) and the angle determined by the central C=C axes 
of the two monomers (bottom) during a molecular dy- 
namics simulation of [TCNE]2~ in dichloromethane. 
Figure 8.: Characteristic frequencies of the TCNE 
monomer and dimer extracted from the velocity auto- 
correlation functions for selected pairs of atoms. 
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FIG. 5: D. Scherlis 
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FIG. 6: D. Scherlis 
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FIG. 7: D. Scherlis 
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FIG. 8: D. Scherlis 



